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A unified energy principle approach is presented for analysing the magnetohydrodynamic (MHD) 
stability of plasmas consisting of multiple ideal and relaxed regions. By choosing an appropriate 
gauge, we show that the plasma displacement satisfies the same Euler-Lagrange equation in ideal 
and relaxed regions, except in the neighbourhood of magnetic surfaces. The difference at singular 
surfaces is analysed in cylindrical geometry: in ideal MHD only Newcomb's [W. A. Newcomb (2006) 
Ann. Phys., 10, 232] small solutions are allowed, whereas in relaxed MHD only the odd-parity large 
solution and even-parity small solution are allowed. A procedure for constructing global multi-region 
solutions in cylindrical geometry is presented. Focussing on the limit where the two interfaces 
approach each other arbitrarily closely, it is shown that the singular-limit problem encountered 
previously [M. J. Hole et al. (2006) J. Plasma Phys., 77, 1167] in multi-region relaxed MHD is 
stabilised if the relaxed-MHD region between the coalescing interfaces is replaced by an ideal-MHD 
region. We then present a stable (fc, pressure) phase space plot, which allows us to determine the 
form a stable pressure and field profile must take in the region between the interfaces. From this 
knowledge, we conclude that there exists a class of single interface plasmas that were found stable 
by Kaiser and Uecker JR. Kaiser et al (2004) Q. Jl Mech. Appl. Math., 57, 1], but are shown to be 
unstable when the interface is resolved. 



I. INTRODUCTION 



The analysis of a toroidally confined plasma begins with finding solutions to the ideal magnetohydrodynamic (MHD) 
equilibrium equation 

VP = JxB, (1) 

where P is the scalar pressure, J = V X B //ip the current density, /xg the permeability of free space and B the magnetic 
field. In configurations that have a continuous symmetry (which we shall call the 2-D case), such as tokamaks and 
reversed-field pinches, or two continuous symmetries (the 1-D case), such as large-aspect-ratio tokamaks and reversed- 
field pinches in the cylindrical approximation, it is possible to find mathematically rigorous solutions to this equation 
by assuming the existence of magnetic fiux surfaces and using the Grad-Shafranov equation. [[1] However, for the 
general case of 3-D equilibria, such as stellarators, which have no ignorable spatial coordinate, magnetic flux surfaces 
are not guaranteed to exist. Indeed, the most generic situation for field lines in 3-D is that they are chaotic and come 
arbitrarily close to any point in a volume. [[5^] 

The problem posed by chaotic fields for the existence of 3-D configurations was discussed by Grad who noted 
in 1967 [[3]] that flux surfaces are only automatic when an axial symmetry is present. However, the Kolmogorov- 
Arnold-Moser (KAM) theorem [[?]] shows that flux surfaces with sufficiently irrational rotational transform, referred 
to as KAM surfaces, may survive symmetry breaking perturbations of a given magnetic field. The situation in a self- 
consistent MHD equilibrium is more complex, but it is reasonable to conjecture that there exist equilibria consisting 
of regions in which the magnetic field is partially chaotic separated by perfect flux surfaces called KAM barriers or 
interfaces. [ 5 ] 

Though regions between the KAM barriers will not in reality be uniformly chaotic, the simplest approach rigorously 
consistent with both ideal MHD and the existence of embedded chaotic regions is to adopt a "worst case scenario" 
and treat these regions as if they were completely chaotic. As chaotic fields cannot support pressure gradients, this 
leads to the stepped-pressure-profile model [[61 [7j], where the pressure and rotational transform steps are positioned 
at fiux surface interfaces and the pressure gradient is zero elsewhere. 

In regions of constant pressure, often called force-free regions in astrophysics, the magnetic field is described by 
VxB = Q:(r)B, where Q(r) is an arbitrary scalar function. By taking the divergence of both sides of this equation 
and using the vector identities V-VxB = and V'[Q!(r)B] = Q;(r)V-B + B'Va(r), we see that = B'VQ;(r) and 
so the function a{r) must be constant along a field line. Therefore, in chaotic regions, the magnetic field satisfies the 
equation 



VxB = aB , 



(2) 
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with a constant throughout the region, such a field being known as a Beltrami field (or linear force-free field). Thus, 
in the stepped-pressure-profile model the magnetic fields in the constant-pressure regions are assumed to be Beltrami, 
but with different a in each region. 

The Beltrami equation is the Euler-Lagrange equation for the relaxed-MHD variational principle of Taylor [[S^] 
for equilibria of plasmas (such as reversed-field pinches) that undergo a strongly turbulent phase during which all 
ideal-MHD invariants except for the magnetic helicity and toroidal and poloidal magnetic flux are broken. Thus 
we shall use the same mathematical framework, the relaxed-MHD energy principle, for modelling both plasmas with 
magnetic-field-line chaos due to 3-D equilibrium effects and turbulent plasma experiments as considered by Taylor. 
Our results in this paper will in fact be more relevant to the 1-D and 2-D systems considered by Taylor as we shall 
mainly work in the cylindrical approximation. The discussion of field-line chaos in 3-D systems above is mainly to 
motivate our stepped-pressure-profile multi-region relaxation model, but recent work [[9j] suggests that multi-region 
relaxation is relevant to the 2-D reversed-field pinches as well. 

As a first step towards solving for the field of these plasmas, Hole et al. [[3 [7]] considered a periodic cylindrical 
model. In this geometry, nested flux surfaces exist everywhere and analytic solutions for the magnetic fleld exist when 
the pressure is constant. Such plasmas form the foundation of this work. 

Hole et al. used a variational treatment to construct equilibrium solutions and analyse their stability in cylindrical 
geometry. [[7|] Using this analysis, they reproduced the earlier conclusion of Kaiser and Uecker [[10]], that stable 
plasmas exist in the case of a single-interface configuration providing there is a jump in the rotational transform, t, at 
the plasma vacuum interface. However, their work also raised questions. First, as part of their stability calculations, 
they compared the stability of a single-interface pressureless plasma with a jump in t to a two-interface plasma with 
the same net jump in the limit that the two interfaces became arbitrarily close. In doing so they found that the 
stability was different in the two cases: the two-interface plasma was unstable when the single-interface plasma was 
not. Second, the current density must be infinite at any interface where there is a jump in pressure. Although this is 
consistent with MHD, it is unphysical in a real plasma, and raises the question of whether the stability conclusions 
hold if interface is replaced by a small region over which the pressure and field vary continuously. 

In this paper, it will be argued that this singular limit paradox is due to the assumption of a relaxed-MHD region 
between the coalescing interfaces, which allows a tearing instability to occur. Instead, the single-interface configuration 
with jump in rotational transform should be compared to a plasma containing a thin but finite ideal region, with 
constraints on the helicity of every field line. If the pressure and magnetic field vary continuously between the values 
on each side of the interface, the infinite current densities of the stepped pressure model are also removed. 

A larger motivation for our study is to remove the singularity in the current density, thereby illuminating some 
limitations in the stability conclusions of Kaiser and Uecker [[10]] ■ We will show that if the interface is resolved as 
an ideal region of nonzero width, the rotational transform profile, t{r), may pass through ideal resonances. If so, the 
stability analysis must consider the structure of the pressure profile in the barrier in the limit it is spatially resolved. 
In some cases, this analysis will force us to conclude that there is no barrier field and pressure configuration that 
is both stable and physically plausible, even when Kaiser and Uecker's analysis showed the plasma could support 
pressure jumps at the interface. 

In order to achieve these goals, it is first necessary to develop a stability analysis that can calculate the stability of 
a plasma consisting of both Taylor-relaxed and ideal-MHD regions. Section ([iT]), derives a stability test for multiple- 
interface ideal-MHD plasmas, based on the calculus of variations and section (III I applies the test to cylindrical 
plasmas. Section (IV I describes how the analysis may be modified to incorporate plasmas with Taylor-relaxed regions. 
Section (|v]) then describes the illustrative two-interface plasma that will be investigated, and the numerical method 
used to generate stable solutions and determine their stability. Finally, section (VI I demonstrates that when the 
resistance of the narrow plasma region is taken to be zero, there is no discrepancy between the one- and two-interface 
results. It will also be shown that the stability of a configuration is related to the form of the pressure profile, not 
just the amount by which it varies between the interfaces. 



II. IDEAL- AND RELAXED-MHD MULTI-INTERFACE PLASMAS 



Our MHD model of an A^-interface plasma draws heavily on two existing approaches. Recently, Hole et al. [[5]], 
developing earlier work by Spies et al. [[ilj] and Kaiser and Uecker [[10 ], determined the stability of an A^-interface 
plasma consisting entirely of Taylor-relaxed regions by constructing an A^ x A^ matrix and calculating its eigenvalues. 
In 1960, Newcomb [[T^] developed a method to analyse the stability of an ideal cylindrical plasma. This section will 
present an extension of the derivation used by Hole et al, using, and comparing with, results of Newcomb where 
appropriate. 

Like Hole et al. , our system comprises A^ nested plasma regions surrounded by a vacuum and enclosed by a perfectly 
conducting wall, as shown in figure ([l]). The plasma regions are denoted by Vi, with Vi at the core, and Vn next to 



FIG. 1: Diagram of plasma regions. 



the vacuum region V. Each plasma region Vi is bounded on the outer edge by the interface Ti and on the inner edge 
by Xi_i, while the perfectly conducting wall, W encases the vacuum. The interfaces are located at r — ri, the wall 
at r = and the width of the plasma is normalised so that rjv = 1. In the following, we adopt the convention that 
when equations refer to a particular region or interface, the region or interface will be given before the equation. 

We employ a variational approach to study equihbrium and stability. The potential energy of a iV-region plasma, 
where each region is described by either ideal or relaxed MHD, can be written 

where P and B are the pressure and magnetic fields and 7 is the ratio of specific heats. We expand W to second order 
in perturbations from its equilibrium value: W — Wq + SW + where 6W, the first variation, is linear in the 

perturbation and S'^W, the second variation, is quadratic (note that we have used the notation of Spies et al. 
where 6^W is twice the quadratic energy functional SW used by [i_13j] or W used by Newcomb The system 

is in equilibrium if the first variation is equal to zero for every possible displacement of the plasma and stable if the 
second variation is positive for every displacement. 



A. Ideal Regions 

To calculate the variation of this integral in plasma regions described by ideal MHD, we introduce the displacement 
vector ^, defined as the displacement of each fluid element as a function of its original position rg. Eulerian variations 
in the pressure, p = SP, and magnetic field, b = (5B, can both be related to ^, provided it is assumed that the plasma 
resistivity and thermal conductivity are zero, [[131 E]] yielding 

P = -(I-VP + 7PV-0, (4) 
b = Vx(^xB). (5) 

Equation ([s]) incorporates a continuum of "frozen-in flux" constraints [[TS]] holonomically. These constraints do not 
apply in the vacuum, so in V we use the representation b = V X a, where is a perturbed vector potential, related to 
^ only at the interfaces through the boundary condition [[13]] 

nXa=-n-^B, (6) 

where B is the unperturbed vacuum magnetic field at the interface. 

Using these expressions, and simplifying the resulting expression for the variation of a volume integral using vector 
identities, we find that 



5W ^ Ylj |-[VP-JxB] dT^ + j a-3 dr^ 
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N 

^=l •'^^ 



P 



(7) 



where n is the unit outward normal at an interface, J — VxB//io is the current density, and \X\ in indicates the 
change in the quantity X across the interface. As the condition 5W = must be satisfied for all possible ^ and a, 
including those that approach a delta function, it is equivalent to 



V, : VP- JxB = 0, 
V : J = 0, 

2^0 



P 



0, 



(8) 
(9) 

(10) 



Since the second variation will only be used to determine the stability of equilibrium solutions, it will be assumed 
that equations (|8]), ^ and (10 1 are satisfied. With these assumptions, the second variation is [|13j] 



N 



b-b 



- ^ X J-b + (V.O(^-VP) + ^P{V-^y' 



dr' 



N 

Jv Mo T^i-Ji^ 



(„.V)(P+^) 



da^ 



(11) 



If is positive for all physical perturbations, the configuration is stable. Finding the global minimum of is 
a non-trivial problem. Instead, as is done elsewhere [[HIII^]] a vector equation will be found which is satisfied by any 
perturbation for which 5^W is a local minimum or maximum. 



B. Relaxed and Vacuum Regions 

As relaxed-MHD regions are by their nature strongly mixed, fiuid displacements are ill-defined in these regions. 
Nevertheless, in the following we discuss how to formulate an energy principle in these regions that is as close as 
possible to that used in ideal-MHD regions. As we can regard the vacuum field as a Beltrami field for which a = 0, 
the discussion applies to the vacuum region as well, where a fortiori ^ cannot be thought of as a fiuid displacement. 

The fundamental representation of the perturbed magnetic field in these regions is b — Vxa with boundary 
condition Eq. As Eq. ([Sl makes it clear that we can take a = ^xB throughout an ideal-MHD region, a unified 

"bed magnetic field is provided by using a gauge such that 



representation for the pertur 



a = ^xB, (12) 



which is consistent with the boundary condition Eq. ^ at the edge of a relaxed or vacuum region and effectively 
defines the projection of ^ perpendicular to B throughout the plasma and vacuum regions (the parallel component 
being arbitrary). This gauge was introduced by Newcomb [|12j] for vacuum fields and we shall henceforth refer to it 
as the Newcomb gauge. 

Newcomb noted that ^ can be singular in this gauge. To see this, let s be a coordinate (e.g. r in cylindrical geometry) 
labelling magnetic surfaces. Then b-Vs = B'V(^'Vs). This can be inverted to give ^-Vs = (B'V)~"'^(b-Vs) except 
at magnetic surfaces where B-V is singular (within the function space of trial functions we allow). If there is such 
a singular surface within the region under consideration, then ^ will diverge in the neighbourhood of the singular 
surface. This will be made more explicit when we consider the cylindrical case in Sec. |IV[ 

Provided we treat these singularities appropriately, the Newcomb gauge will allow us to unify the discussion of 
ideal, relaxed and vacuum regions. However, before discussing the application of the Newcomb gauge to relaxed MHD 
we first review the standard approach of using an arbitrary gauge. 

The multiregion [|S] |T5]] generalization of the relaxed-MHD energy principle is to minimize the total potential 
energy W Eq. ^ under the constraints of constant magnetic helicity 



K, = ^ [ A-B (^r^ (13) 

'Vi 



2fiQ 



and constant PiV^ in each relaxed region, of volume Vi, where P = P^ is constant throughout the region (i.e. VP = 
in Vi). Note that we have defined K with a prefactor of 1/2/io, so the helicity as commonly defined [[S^] is 2fioK. 
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The pressure constraints can be implemented holonomically [fl6]], while the helicity constraints are implemented by 
introducing Lagrange multipliers ai, and replacing the constrained problem SW = with that of finding an extremum 
of a magnetic "free energy" [[17]], 



N 



i=l 

subject only to the constraint that variations a not vanishing on the boundary satisfy Eq. ([6|. 
The first variation is [ [TTIITB] ] 

N+l 



(14) 



SF = y2 a.(VxB - a^B) dr 

i=i 

N 



2mo 



Setting 8F — for all possible a and ^, we have the equilibrium equations 



V, : VxB = a,B, 
52 



P 



2^0 



0, 



(15) 
(16) 

(17) 
(18) 



Eq. (17) being a Beltrami equation Eq. ([2| in each relaxed region, with oat+i = 0. 

It can also be shown for Beltrami fields (cf. [[H]] and provided the trial function space of ^ is such that 6oLi — 
and SVi — [[16 ]) that the constrained second variation S^W is the same as the unconstrained second variation S^F 



S^F = [ — (b-b- a,a-b) dr 

A' 



(n-V) 



2/io 



da-" 



(19) 



We now compare the relaxed-MHD and ideal-MHD variations, S^F and 6"W respectively, to show they are equiv- 
alent in the Newcomb gauge Eq. (12 1 when evaluated on a relaxed equilibrium, i.e. when the Beltrami equations 
Eq. (17 1 are satisfied. With this assumption the ideal first variation Eq. ^ implies VP = and the interface force 
balance equations Eq. (18 1, as expected. The assumption 6Vi = allows us to set V-^ = (see Sec. Ill I. Then, 
recognizing that ^xJ-b = aia-h/^o, we see that Eq. (11) is indeed the same as Eq. (19) for a Beltrami equilibrium. 



Because of this equivalence we shall simply take it as understood that S W means S'^F in a relaxed-MHD region. 



C. Kinetic Energy Normalisations 

When searching for local minima of S^W, two problems are immediately encountered. 

First, since S'^W scales with the square of the amplitude of ^, a finite minimum need not exist (and, if it does, it is 
the trivial minimum S'^W = attained when ^ = 0). To compensate, we normalise S^W by assigning a nonzero value 
to a positive-definite quadratic functional, M[^]. For the purpose of determining stability from the energy principle, 
the choice of normalisation is not unique; apart from being quadratic and positive it is necessary only that M be 
finite for all admissible functions ^. In this section we discuss criteria for choosing M, interpret it in terms of a kinetic 
energy, and make a simple choice for use in the following sections. 

The second problem is that the values of A = S'^W/M may approach a lower bound that is not attained by any 
physical perturbation. The class of perturbations considered will therefore be extended to a generalised function space 
that includes any vector function with the following properties: 

1. The function may be approximated to any degree of accuracy by a physical perturbation. 

2. The value of S^W for the vector function is finite and approached by the values of 6^W for approximating 
perturbations. 
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If is positive for all ^ satisfying these requirements, then it must be positive for all physical perturbations, and 
so the configuration is stable. If a non-physical perturbation exists for which S'^W is negative, the second requirement 
dictates that there also exists a physical perturbation with negative (5^W.[[T2]] 

To find the minimum value of S^W, subject to the constraint of fixed M, the technique of Lagrange multipliers is 
used, setting SL = 0, where the Lagrangian L = — vM .[\1V] As L is a symmetric quadratic form with respect 
to ^ [[18 ], if 5L = 0, then i = 0. Therefore, when 5L = 0, the value of i' is the same as the value of A = S'^W/M. 
The stability of the configuration can therefore be determined from the sign of the Lagrange multiplier, which will 
hereafter be referred to as A. 



A positive eigenvalue can be interpreted as the square of the frequency of a normal mode, A 



while a negative 



eigenvalue, which implies instability, is the square of the growth rate, A 



-7 



In this interpretation 5 AM is a model 



kinetic energy. As it chosen for mathematical convenience, [j/7 is not necessarily an accurate estimate of a physical 
frequency /growthrate . 

Several choices for the model kinetic energy are reviewed in [[H]]. In ideal MHD we can in fact normalise with the 
physical kinetic energy factor M = //? |^ P dr^ where p is the mass density. However, to remain closer to Newcomb's 
[[12^] stability analysis it is convenient to replace p with an anisotropic tensor density p = pnn, so that, in an 
ideal-MHD region Vi, we use 



p{n.if dr^ 



(20) 



which has been used for instance in the PEST2 code [[201 EI]]- This defines not only a normalisation but a norm 
for the function space in which admissible solutions in Newcomb's formulation exist. In this space the linearised 
force operator is self adjoint and the integration by parts required to obtain Newcomb's ideal-MHD 1-D scalar energy 
principle, Eq. (15) of [^], or Bineau's 2-D generahsation of it [UHl HH HH US]], can be carried out. 

In the case of relaxed MHD and the vacuum, dt£, cannot be interpreted as a physical velocity and Mn is not a 
physically natural norm. In fact, as it is crucial that the helicity be finite in relaxed MHD the natural norm is one 
based on the helicity. Thus we define the helicity norm for use in a relaxed-MHD or vacuum region 



MnM = I / /i|xB-Vx(^xB) dr^ 



(21) 



where h is an arbitrary factor. 

Choosing h so MH,i has the same dimensions as Mn,^ we can define the total normalising factor M as a sum over the 
regions, using A/n.i or Mjj^i as appropriate. However, provided the plasma is stable against perturbations that do not 
move the interfaces, we may effect further simplification by taking p and h to be Dirac delta functions concentrated 
on the interfaces and use 



JV 



(22) 



(If the above assumption is violated then this normalisation will give A = —00.) 

In the following we use Eq. (22) for normalisation and, with p — \ and h= \, use Mn.; and A/h.; as norms to select 
the appropriate function space for ^ in each region: The "Newcomb" norm, MN,i, restricts the trial function space for 
n-^ to squarc-integrable generalised functions, while the helicity norm, Afu.i, allows more divergent functions. This 
will be made more explicit in Sec. |III[ 

In evaluating 5L, we consider the equilibrium quantities as fixed, varying ^ only. Using vector identities and 
equation (11), we find 



5L = 



I ^i- {Bxj - Jxb - V{i-VP) - V[7P(V.O]} dT^ 
[[B.b+(^.V)(B2M1 ^ 



(5^-Bxj dr^ 



N 

E 

1=1 



Mo 



- A(^-n) da^ 



(23) 



where j = Vxb//io. Since 6L = for any possible perturbation, including those that approach delta functions, we 
may conclude that 



V^: Bxj - Jxb - V(^-VP) - V[7P(V-a] = 0, 



V 



Vxb = 



2/io 
0. 



Bb 



Mo 



- A(^.n) = 



(24) 
(25) 
(26) 
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Equation (24 1 is a system of three second order partial differential equations, and therefore does not have a unique 



solution. In order to proceed further with this analysis, it is necessary to specify a geometry for the configuration. 



III. CYLINDRICAL IDEAL-MHD MULTI-INTERFACE PLASMAS 
A. Solutions in Case of no Singular Surfaces 

This section will apply the stability analysis to periodic ID cylindrical plasmas with coordinates (r,9,z), in which 
the magnetic field and pressure are functions of only the radial coordinate, r. By imposing a periodicity constraint of 
length L in the z direction, the model is a simple analogue of a large-aspect-ratio, toroidally symmetric plasma with 
circular flux surfaces. In this analogue, 6 and 2ttz/L correspond to the poloidal and toroidal angles, and L — 2ttR, 
with R the major radius of the toroidal plasma. 

Like Hole et al. [[Bj], we next decompose ^ into a Fourier series with components: 



$.m,k{r, 9, z) = [Cm,fc,r(?')f -|- £,mM,e{r)0 + £.7nM,z{r)2] cxp[i{m9 + kz)]. 



(27) 



Here, the first two subscripts refer to the wave number in the 9 and z directions respectively, the third subscript 
refers to the geometric component, f , 6 and z are unit vectors, k = —2'!m/L, and m and n are integers. The second 
variation is decomposed as a sum of uncoupled terms 



(28) 



m.k 



where S^Wm.k is the second variation evaluated for the m, k component. Thus we can consider each Fourier component 



individually, restricting the trial function space to a single value of m and fc, so that Eq. (24 1 reduces to a second 
order linear ODE. 



Under the Fourier decomposition of equation ( 27 1 , the dot product of equation ( 24 1 with B reduces to 



B.V[7F(V.e„,fc)] -0, 



(29) 



where we have used B-VP = and equation (l8|. The quantity 7-P(V'^m,fe) is therefore constant along field lines, 
and by extension on magnetic surfaces in the cylinder. As the pressure is constant along magnetic surfaces, V'^m,fc 
must also be constant. However, V-^m.fc is a function of r multiplied by exp[i(m0 -I- kz)] and so if the divergence is 
to be constant on a magnetic surface, ^ •$,m,k = 0; except \i m = k = Q. Therefore, ignoring the case m = fc = 0, the 
normal modes $,m,k must satisfy 



B X ]rn,k — J X h^^k 
^•inr,k = 0, 



(€m,fc-V) 



2mo 



Bb 



Mo 



A(^m,fe(r)-n) = 0, 



(30) 
(31) 

(32) 



where hm,k = Vx(^,„,fcXB) and j„i,fe = Vxb„,.fe. 

First we solve for ^m,fc in the plasma. Using the 9 component of equation (30 1, along with the requirement $,m^k is 
divergence free, we are able to eliminate S,m,k,e{r) and im,k.z{r), yielding 



dr \ dr 



.(r) 



(33) 



where 



f{r) 
9{r) 



[mBe{r) - krBz{r)Y 



r{k 



[mBe{r) + krBz{r)]'^ 2Be{r) d 
r r dr 



[rBe{r)] 



d 

dr 



m^Beirf - Pr'^B^ir)'^ 



(34) 



(35) 
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This is equation (23) of Newcomb [fT^], which was derived by first eliminating ^m,k.e{r), fm,fe,z('') and then taking 
the variation, rather than the opposite order presented here. Newcomb uses equation (l8| to rewrite g{r) as 



dP [niBeir) + krB^{r)]'^{Pr'^ + - 1) 



^2^2 _|_ ^2 ^j. 



Equation ( 33 1 has a singular point at r = and whenever 

krBz{r) + mBg(r) 



0, 



0, 



or equivalently 

k i(r) 

h 

m R 

where t(r) = i(r)/(27r) is the rotational transform, defined in cylindrical geometry as 

RBe{r) 
rB,{r) ■ 



*(r) 



(36) 



(37) 



(38) 



(39) 



Surfaces corresponding to values of r at which equation ( 37 1 is satisfied are also known as resonant surfaces since on 
these surfaces the phase of the displacement is constant along magnetic field lines, i.e. k-B = 0, where k is the wave 
vector. It is at these surfaces that the operator B-V discussed in Sec. IIB| is singular in the trial function space m, k. 
Initially it will be assumed that there are no solutions to equation (37 1 in a given region. The case when singular 
points are present will be treated separately. 

Next, we choose a representation for ^m,k on an interface X^. To evaluate the interface terms, hm,k is expanded in 
equation (32) and the expressions for ^z{f) and Ce('r) found when solving equation (24) are substituted in, yielding 

[[k^r [Bzir^ + Be{r)^] ^r)]], [/(r)e/(r)] 



(40) 



Since equation (33) has two linearly independent solutions, it is necessary to find boundary conditions that allow easy 
evaluation of equation (40), ultimately by converting the N equations to a matrix. Let Xi be the value of ^m,fc.r (^) 
on the ith interface, which is located at r = r^. Ignoring for the moment the first plasma region and the vacuum, the 
solution to equation ( 33 1 that satisfies the boundary conditions 

Cm,fc,r(^i— l) — 

may be written as 



1 <r < Ti, 

where Ci.i'" and Ci,2(^) are two independent basis solutions that satisfy the boundary conditions 

C,,i(r,„i) = l C„i(r-z) = 0, 
0,2(^.-1) =0 0,2(r,) = l- 



(41) 
(42) 

(43) 



Substituting Eq. ( 42 1 and its derivative into equation ( 40 ) , we see that the interface equations depend linearly on the 
Xi. In the case of the outermost interface, the requirement on the wall W : n-b = provides the boundary condition 
^m,k,r{'rw) — and the surviving basis solution Qn+i is normalised to 1 at rjy. In the first plasma region, the single 
basis solution is chosen so that it remains finite at the singular point r = and is normalised to 1 at ri . This solution 
satisfies the boundary conditions [[T^] 



Ci(o) : 
Ci'(o) 



0, |to|^1 

: 0, ItoI = 1 



(44) 



Using the solutions and Cj_2, the N interface equations given by equation (40) may be written in matrix form 

AX = AX (45) 

where 

X = {Xi,...,X^}. (46) 

The stability of the m, k component of ^ can then be confirmed by checking that all eigenvalues of the N x N matrix 
A are positive. 
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B. Ideal Solutions with Singular Surfaces 



Up to this point, it has been assumed there are no resonant surfaces in the plasma region Vi- Under this assumption 
the solutions to equation (33 1 will be finite and continuous, corresponding to physical perturbations. However, if there 



is a resonant surface in the region, the solutions are not guaranteed to be finite, and so it is necessary to check that, in 
an ideal-MHD region, the solutions meet the two conditions for test functions given in section ([iT]): the function may 
be approximated to any degree of accuracy by a physical perturbation and the value of 5^W for the vector function 
is finite and approached by the values of for approximating perturbations. 

Let us assume there exists one singular point in Vi , that satisfies equation ( 37 1 . By taking a series expansion of 
the functions f{r) and gir) in equation (33) about subject to the condition that equation (37) is satisfied, it can 
be shown that, to leading order in r — r^, /(r) = a{r — and g{r) — j3 [pL2j], where 



Making the change of variable x = r — Vg, equation (|33|) may then be approximated by 



d ( 2d^r{x) 

a— X — 



dx 



dx 



The solutions to this equation have the functional form \x\^ , where 



V = 



4± 



(47) 
(48) 

(49) 
(50) 



These solutions form the leading terms of the full Frobenius expansion of about the singular point. If the two rj 
roots differ by an integer, the full expansion may include logarithmic terms. [[24]] This is particularly relevant to this 
work, because when P'{r) = the two solutions to equation (50) are and —1. It turns out that in the special case 



of a Beltrami field the logarithmic terms vanish, however this is not the case for constant pressure plasmas in general. 

If a + 4/3 < 0, the two values of rj are complex, and so the solution oscillates. It is also the case that the frequency 
of oscillation increases as the singular point is approached. Such solutions are associated with instability [[I2]]) [[25]] 
and so a first requirement for stability is that a + 4/3 > 0. Using the definition of * and the resonance condition, the 
inequality may be reduced to Suydam's condition 



rB,(r)2 /dlog(*) 



8i?2 



dr 



dP 
dr 



> 0. 



(51) 



If Suydam's condition is satisfied, the two values of rj are real. One of the solutions, known as the small solution, 
rjs > —1/2, diverges slower than l/(r — Ts)^^^ as r — > r^, the other, known as the large solution, rji < —1/2, diverging 
faster than 1/(7- — rsY^^. Because the singular point at rg is excluded, these solutions are not ordinary point wise 
solutions of Newcomb's Euler-Lagrange equation, but are better regarded as generalised functions [[26']. From this 
point of view there are actually four solutions in the neighbourhood of r^: two small solutions and two large solutions, 
depending on whether they have support on one side or the other of r^ (or, equivalently, small and large solutions that 
are odd and even about to leading order). However, the two large solutions are excluded in ideal MHD because 
the norm A/n,! defined in Eq. (20) is infinite. Thus, in ideal MHD we still have only two independent solutions as in 
the non-singular case, but they are discontinuous at r^, as noted by Newcomb [[H]]. 



Using a similar construction to that used when no singular points are present, two basis solutions to equation (33) 
are defined according to the boundary conditions 



Ci,i,s{ri-i) = 1, lim^^r- CiA,sir) = a\r - r^C^, 0,1, s(f) = for r G {rs,ri], 
0,2,s(''i) = 1, l™r_»r.+ Ci,2,s('') = b\r - r^C^, Ci,2,s('") = for r e [r^-ijr^), 



(52) 



where r]s is the value of 77 corresponding to the small solution and a and b are constants determined by the interface 
boundary condition. As before, the subscript i denotes the plasma region, and the second subscript distinguishes 
the two independent solutions. The third subscript indicates the presence of resonant surfaces. Equation (42) still 



applies, but with and replaced by Ci,i,S and 0,2,5', respectively. In contrast to the case where there are no 
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resonant surfaces, where (/{ri) is a function of Xi-i, Xi and Xi+i, Q s'(''i) depends only on Xi. As a consequence, 
two adjacent off-diagonal entries in the tri-diagonal matrix A are removed, allowing it to be separated into two smaller 
matrices for which the eigenvalues may be calculated separately. 

Finally, if a singular point exists in the inner plasma region, the boundary condition that Ci{r) is finite at r = is 
replaced by the requirement that the solution is small as r — > Ts from above and vanishes in the interval < r < r^. 
This modification alters the eigenvalues of the matrix A, but not its structure. 



IV. TAYLOR-RELAXED AND IDEAL-MHD MULTI-INTERFACE PLASMAS IN CYLINDRICAL 

GEOMETRY 



This section will describe how to adapt the analysis to plasmas consisting of both ideal regions, in which the 
analysis of the previous section is applied, and Taylor-relaxed regions, where the analysis of Hole et al. is used. 
The analysis of Hole et al. may only be applied to Beltrami fields, where VP = and VxB = aiB. The stability 
analysis developed in this paper makes no assumption on the form of the equilibria in ideal-MHD regions, which can 
be either Beltrami states or regions with finite pressure gradients. 



Suppose the equilibrium field in Vi is Beltrami. Equation (24 1 reduces to 



^ ^ Bx[Vxb„,fc] _ a,B X hra,k ^ Q ^gg^ 
* ' Mo Mo 

where is the Lagrange multiplier of the Beltrami field. Since this equation does not involve the perturbation ^, it 
is possible to solve for hm.k directly 

T'., : Vxb^^fc = aib™,fc ttm^feB, (54) 
where a^.k has the same 9 and z dependence as b„j fc, i.e. exp[i(m0 + kz)]. Its r dependence is determined by 



taking the divergence of both sides of Eq. (54l, giving B-Vam^fc = i{mBg/r + kBz)am,k = 0. Thus am.fc vanishes 
except in the neighbourhood of a singular surface, where it can have the form of a Dirac delta function with arbitrary 
amplitude, corresponding to a current sheet at the singular surface. 



Apart from the current sheet term, Eq. (54) is the same as the relaxed-MHD volume equation used by Hole et al. 



Put another way, am,k = in relaxed MHD. Thus, using the normalisation Eq. (22), the existence of current sheets 
at resonant surfaces is the only feature that distinguishes linearised relaxed-MHD perturbations from ideal-MHD 
perturbations on a Beltrami equilibrium region (as expected from our general discussion in Sec. [H]). It is the current 
sheet that allows the small solutions on either side of a singular surface to be disconnected in ideal MHD , screening one 
side from the other. In ideal MHD, field-line reconnection is forbidden by the frozen-in flux condition, so the current 
sheet must form to prevent the tearing-mode island (of arbitrarily small amplitude in the linearised approximation) 
that forms in relaxed MHD [[8 ]. 

In order to modify the ideal analysis of Sec. |III| to include Beltrami regions, it is only necessary to change the 
boundary conditions used to define the solutions Q i s (^) 2 s ('')• Although the perturbation function ^m.k,r 



has no physical meaning in the Taylor analysis, we extend the definition for mathematical convenience as in Eq. (12). 
Using the identity b-Vr = B-V(£-Vr), the radial components bm^k,r{f) and £,m,k,r{T) are related by 

^ra.kAr)- ^^^^^^ . (00 j 

As the r-component of b must be continuous in r for it to be possible to satisfy V-b = 0, bmkrir) is an even 



function in the neighbourhood of a singular surface (i.e. asymptotically as r —> r^). Thus Eq. (55) implies that 



for relaxed-MHD perturbations, only the asymptotically-odd large (,m,k.r solution and the asymptotically-even small 



solution are allowed. (This is consistent with the hclicity norm, Eq. (21), defined as a principal part integral, as the 
odd large solution has finite norm but the even large solution does noTT) 

Thus, as in the ideal case discussed in the previous section, there are also only two independent solutions in a 
relaxed-MHD region containing a singular surface. Using these results to continue across singular surfaces, relaxed- 
MHD basis solutions are then defined so the boundary conditions 

Ci,i:s(?'j-i) = 1 Chi^Bin) = 

Ct,2,B{ri-i) = Ci,2,B{ri) = 1, (56) 

are satisfied. The analysis then proceeds as if there were no singular points in this region, with Cj:,i,b(^) and Cj:,2,b(^) 
taking the place of Ci.i(^) and Ci.2{i^) in forming the matrix A. 
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V. TWO-INTERFACE PLASMA 



This section describes the numerical techniques used to apply the analysis to a two-interface plasma. In the 
plasmas considered, the central plasma region is assumed to be a force-free Beltrami equilibrium state while the edge 
region between the internal interface and the plasma-vacuum interface is a pedestal region which can support a finite 
pressure gradient (provided we use ideal MHD in this region). To describe the equilibrium in the edge region, the 
rotational transform and pressure profile were specified and the pressure variation accounted for through variation in 
the magnitude of the magnetic field. For such configurations, the magnetic field takes the form 



V : 



Be{r 

Bg{r 
Bg{r 



= CiJi(Q;ir) 
= CiJo(air) 

= BX/r 
= B, , 



(57) 



(58) 



where ci, the field strength coefficients, B^ and Sj, 
Several factors do not influence stability. First, 



and the Lagrange multiplier, ai, are constants, 
rescaling the magnetic field by a factor, a, does not change 
the ODE for S,r{r) [equation (33l] providing that the pressure profile is also rescaled by a factor of . As such, 

B^ — J {B'XY + {BYY ™3'y be set to a convenient value, in this case JJIq. Second, changing the sign of both ai 



and t(r) changes the sign of Bg(r) everywhere. Since equation (33 1 is unaffected by changing the sign of k and Bg(r) 
together, only i{ri) > need be investigated, with k spanning both positive and negative values. Finally, the value 
of R plays a limited role in stability calculations, as it determines which values of fc = —n/R, n an integer, should be 
considered and elsewhere appears only in the resonance condition, equation (37 1. Without loss of generality, we select 
R — 1 and make k a continuous variable. Our stability tests will therefore check if a configuration with a particular 
i/R profile is stable for all values of R. If the rotational transform profile is continuous, the equilibrium is defined by 
the rotational transform profile in the second plasma region and the pressure profile, which may be used to determine 
the value of ai and the vacuum field components. 

It was proved by Newcomb [[H]] that to demonstrate a configuration is stable, it is sufficient to show that it is 
stable for all values of k when m = 1, and for m = 0, A; ^ 0. We have set m to 1, and scanned k over the range k^nin 
to /cniax in steps of 0.1. As it was observed that minimum values of A typically occur near = 0, /c = — ti or fc = — 12, 
the values /cmin = —20 and fcmax — 20 were chosen. If the configuration was found to be stable, a finer scan, with a 
step size of 0.01 was then performed around the least stable k values. We also set to to and scanned from k = —0.5 
to fc = 0.5 in steps of 0.01. 

Three methods of stability analysis are of interest: Taylor- Taylor (possible only when the field in the edge region 
is Beltrami), Taylor-ideal and ideal-ideal, where the names refer to the type of analysis used in the core and second 
plasma regions respectively. Solutions for the functions Q were found using different methods for the two regions 
i = I and i — 2. 

In the inner plasma region, the analytic solutions for hm^k and equation (551 provided an exact expression, 
two arbitrary constants were chosen so the function met the appropriate boundary conditions: Ci 
(ideal analysis with no singular points and Taylor analysis) or approaches Vg for the small solution (ideal analysis 
with singular point at Ts), and Ci(ji) — 1- 

In the second plasma region, Mathematica[[27 ] was used to find the solutions numerically. The methods used were 
different depending on the analysis used, and whether there were singular points in this region: 



The 

is finite at r = 



(i) Ideal analysis with no singular points: 

boundary condition 



the numerical solution, C2^i(r') to equation (33 1 that satisfies the 



/-N 
^2,1 



(59) 



was found. The value of (^21 '(''2) is arbitrary, and was set to unity for convenience. C2,i('') is defined by 



C2,i(r) 



C2^l(0 



(60) 



Using a similar method, C2,2(?') was also determined. 

(ii) Ideal analysis with singular point at rgi having determined the location of the singular point numerically, 
solutions either side of the singular point must be found. For r < r^, the solution (^21 s(^) to equation (33 1 that 
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FIG. 2: Stability comparison of Taylor- Taylor and ideal-ideal analysis for m = 1 modes. Panel (a) shows the rotational 
transform profile used and panel (b) a dispersion curve for the Taylor- Taylor (dotted) and ideal-ideal (dashed) analysis. 



satisfies the boundary conditions 



Ci',i.sirs-e) = e^^ C?.s'ir.-^) = -Vse^''\ (61) 



was found, where rjs is the small value for 77 given by equation (50 1. The value of e was chosen so that the 
approximation ^(r) « |r — rs\^^ is valid in the region — e < r < r^. To do this, solutions were found for a 
sample configuration for various values of e and the values for ^(ri) compared. It was found that for e < 10~^ 
there was little difference in the values, and so e = 10^^ was used. The basis function C2.i,s{'r) defined in section 
( |III[ ) was then determined by dividing the function by its value at ri. Using a similar method, C2.2,5(?') was 
determined from the function that satisfies the boundary conditions 

C^.2,si^s + e) = e^^ C^,2,s'irs + e) = ^,e^^-\ (62) 

(iii) Taylor analysis: although the analytic solutions could have been used, in this work the solutions were found 
numerically. The boundary conditions satisfied by C2,i,b(?') and C2.2,s('') were converted to boundary conditions 



for the radial component of the magnetic field using equation ( 55 1 . The equation V X B = 0:26 was then solved 



numerically with these boundary conditions. Equation ( 55 1 was used again to convert back to solutions for 
C2,i,B{r) and C2,2,s(f)- 

In all cases, the solutions for Ci('')i C2('") and Cvi^) were then used to calculate the matrix elements of A. The 
eigenvalues were found, and plotted as a function of k. 

To validate the analysis, configurations with Beltrami fields in each plasma region were analysed under both 
the Taylor- Taylor and ideal-ideal treatments. The Taylor- Taylor analysis was checked by comparing an eigenvalue 
dispersion curve with one generated by the treatment used by Hole et al. [[7]]. As the two methods are equivalent 
when no singular point is present, the ideal-ideal curve must agree with the Taylor- Taylor result for these values of k. 
Figure ([2]) (a) shows the rotational transform profile of the benchmark configuration with ri = 0.5, t{ri) = 1.9 and 
*(''2) — —1-4, and figure ^ (b) shows the m — 1 dispersion curve. These dispersion curves show the values of the 
two eigenvalues. A, of the matrix A as a function of k. If one of the eigenvalues is less than zero, the configuration is 
unstable. There is a resonant point in the edge region whenever k > 1.4 or A; < —1.9. As expected, figure ([2| shows 
that the two methods do not agree for this range of k, with the ideal-ideal analysis giving the more stable eigenvalues. 
If there is no resonant surface in the edge region, the two methods agree, with the exception of a small deviation, not 
visible on this scale, corresponding to the presence of a singular point in the core region. 
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FIG. 3: Rotational transform profile of reference configurations used for vanishing interface separation studies. Note the the r 
axis has been distorted to magnify the second plasma region. Panel (a) shows a configuration for which the rotational transform 
passes directly from t(ri) = 1 to t(r2) = 4, and panel (b) shows a configuration for which the rotational transform passed 
through infinity. 



VI. LIMIT OF VANISHING PEDESTAL WIDTH 



To resolve the vanishing interface-separation paradox, the single-interface dispersion curve was compared to the 
dispersion curves produced by the Taylor- Taylor and ideal-ideal treatments of the two-interface plasma in which the 
width of the edge region, Ar ~ ^2 — ri, is small. As observed by Hole et al. [[7]] when the perfectly conducting interface 
is replaced by a small Taylor-relaxed region (the Taylor- Taylor treatment), the two-interface plasma is unstable where 
the single-interface plasma is not. It will be demonstrated here, that when the edge region is assumed to be ideal, the 
one and two-interface results are in agreement. 

We base our reference example on a plasma found stable by the analysis of Kaiser and Uecker: a single interface 
plasma with a Beltrami field, r^, = 1.5 and rotational transform that jumps from 1 to 4 at the plasma- vacuum 
interface (in the parameters of Kaiser and Uecker, /i — 1.435 and S = 0.540).[ lOj] When the interface is replaced 
by a finite edge region, there are two possibilities for the rotational transform profile in the edge region. Either it 
passes from t(ri) = 1 to *(r2) = 4 directly (hereafter referred to as the fundamental), as shown in figure ^ (a), or 
it diverges to infinity (hereafter referred to as the harmonic), as shown on panel (b). In the example used, the field 
in the edge region is Beltrami, but the conclusions hold for any rotational transform profile with these properties. It 
will be shown that in each case it is mathematically possible to create a configuration with finite edge width that is 
stable, however such configurations are physically unrealistic. 

Our initial stability analysis of plasmas with these rotational transform profiles has a pressure profile that is zero 
everywhere. The eigenvalue dispersion curves produced using both the Taylor- Taylor and the ideal-ideal treatments 
for the m = 1 mode are shown in figures Q, for the fundamental, and ([s]), for the harmonic. Since the m — mode 
was stable in the limit fc — > 0, the m = I dispersion curves provide a sufficient condition for stability. In each figure, 
the two panels show the same dispersion curves on different scales. Panel (b) focusses on the region around marginal 
stability, and includes for comparison half the single-interface eigenvalue for a configuration in which t jumps from 
t = 1.9 (inner side of interface) to t = —1.4 (outer side of interface). The factor of 2 arises from the normalisation, 
equation (22 1, since the two-interface configuration contains twice the number of interfaces. The expression for the 
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FIG. 4: Eigenvalue dispersion curve for m = 1 modes of configuration with rotational transform profile show in figure fSjl (a), 
ri = 0.999, Tw = 1.5 and P(r) = 0. Panel (a) shows the Taylor- Taylor (dotted) and ideal-ideal (dashed) eigenvalues. Panel (b) 
shows the same eigenvalues, focussing on the region around marginal stability. Also shown is half the single-interface eigenvalue 
(sohd). 



single-interface eigenvalue, I, found by rearranging equation (40 1 using equation (10), is 



^ ^ h{ri)Ci'{ri) - fv{ri)Cv'{r,) 2fcViAP ^ 

where ri is the location of the interface and AP is the pressure on the inner side of the interface minus the pressure 
on the outer side. 

Looking first at the fundamental, under the Taylor- Taylor analysis, one of the eigenvalues agrees with half the 
single-interface result, while the second eigenvalue is very positive for k < —4 and fc > — 1 and very negative for 
—4 < A: < — 1, the latter result being responsible for the conclusion of instability. The range —4 < k < —I also 
corresponds to the values of k at which there is a resonant surface in the edge region. Outside this range, the two 
methods produce the same eigenvalues, as is expected. Within this range, neither ideal-ideal eigenvalue agrees with 
the single-interface case, and one eigenvalue is negative for —4 < fc < —1.5, meaning the configuration is unstable. 

For the harmonic, we find that one Taylor- Taylor eigenvalue agrees with half the single interface result, while the 
other is very positive for —4 < fc < — 1 and very negative for fc < —4 and fc > — 1. The range of fc values for which 
the Taylor- Taylor analysis is negative again corresponds to the fc values at which there is a resonant surface in the 
edge region. The ideal-ideal analysis has negative eigenvalues for — 1 < fc < 0. 

To understand why the ideal-ideal analysis may produce negative eigenvalues when the single-interface result is 
stable, the expressions for the the eigenvalues, at fc values with resonant surface in the edge region in the limit 
Ar 0, 

, /i(ri)Ci'(ri) 2fcV| 

2fc2r| fv{r2)Cv'{r2) 
= fc2.2 + ^2 iPi-^) - Pirs)] , (65) 
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FIG. 5: Eigenvalue dispersion curve for m = 1 modes of configuration witli rotational transform profile show in figure (jsjl (b), 
ri = 0.999, = 1.5 and P{r) = 0. Panel (a) shows the Taylor- Taylor (dotted) and ideal-ideal (dashed) eigenvalues. Panel (b) 
shows the same eigenvalues, focussing on the region around marginal stability. Also shown is half the single-interface eigenvalue 
(sohd). 



will be examined. The expressions are equations (93a) and (93b) of Newcomb [[12 ], converted to the notation used 
here, and multiplied by l/(A:^r| -I- m^). If the configuration is stable, both eigenvalues must be positive, and the 
expression may be rearranged to yield 



'2 



2 



, > P(r,(fc)) > P(ri) - ' " 'on , (66) 



where it was noted that P{r2) — 0. The pressure at the singular points is written in the form P{rs{k)) to emphasise 
that the pressure bounds are determined not by the position of the singular point, r^, but by the value of k with 
which it is resonant. Put another way, if we know that the rotational transform varies between the values ti and 



and which harmonic it follows, equation (66 1 gives information about the form a stable pressure profile must take, 
without it being necessary to know the precise location of any of the singular points. 

To illustrate the significance of these requirements, figure ^ (a) shows the lower and upper bounds of equation ( |66| 
as a function of k when ti — 1, i2 — ^ and P{ri) — 0, the configuration used to produce the plots of figure ([3]). This 
figure reveals that if the rotational transform is resonant with k values in the range —4 < fc < —1.5, the plasma may 
only be stable if the pressure is negative at the radial position at which these resonances occur. As the fundamental 
must be resonant with this range of k values in the edge region, this class of configurations may only be stabilised by 
allowing negative pressure. Since the configuration analysed in figure (|4| had zero pressure everywhere, this explains 
why the configuration was unstable over this range of fc. Similarly, the pressure must be positive wherever resonances 
with fc values in the range —1 < fc < .5 occur. Since harmonic configurations must have such resonances, the harmonic 
may only be stable if the pressure in the edge region reaches a greater value than that on the interface. 

When the stability analysis was repeated for the pressure profiles shown in figure ([3]), chosen as example profiles 
that satisfy the pressure bound, the configurations were found to be stable for m = 1 modes and m = in the limit 
fc 0. It is therefore mathematically possible to construct a pressure profile that stabilizes the configuration with 
*i = 1, *2 = 4 and = 1.5, in agreement with the conclusion of stability reached by Kaiser and Uecker. However, to 
achieve this it was necessary to either allow the plasma to have negative pressure, or a very narrow pressure maximum 
at the edge, both of which are physically unlikely. 
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FIG. 6: (a) Stable region (shaded) of P{rs{k)) under ideal-ideal analysis for configurations with ti — 1, t2 ~ 4, ro = 0.999 and 
= 1.5. (b) Stabilising pressure profile for fundamental, shown in figure ([3| (a), (dashed) and harmonic, shown in figure ([3| 
(b), (dotted) configurations. Note that the r axis has been distorted to magnify the second plasma region. 



Rearranging inequality (66 1 reproduces the inequality / > 0, the single-interface stability condition of equation 



(63 1. Therefore, under the ideal-MHD treatment, the stability of a two-interface plasma in the limit Ar may be 
reconciled with the stability of the corresponding single-interface plasma. Specifically, if the single-interface eigenvalue 
is greater than zero for all values of k, then: (i) if there is no resonant surface between the two interfaces, one of the 
eigenvalues will agree with half the single-interface eigenvalue, while the other is extremely positive, and (ii) if there is 
a singular point between the two surfaces, it is possible to construct a pressure profile using inequality (66) for which 
both eigenvalues are positive. If the single-interface configuration is unstable, then the two-interface configuration 
will be unstable under both Taylor- Taylor and ideal-ideal treatments. 

Although we have shown that it is mathematically possible to reconcile the one and two interface stability results, 
the requirements placed on the pressure of the two interface plasma may be physically implausible. It is therefore not 
sufficient to simply apply the single interface analysis of Kaiser and Uecker. It is also not necessary to have a detailed 
knowledge of the plasma in the edge region. Rather, one must consider the field and pressure profiles of the core and 
vacuum regions and the resonances in the edge region, and determine if the restrictions they place on stable pressure 
profiles are physically acceptable. 

Although it was not possible to construct a physically reasonable configuration from the single interface plasma 
considered here, this is not true for all single interface plasmas considered by Kaiser and Uecker. Stable reverse field 
pinch configurations, with rotational transform profile similar to that of figure ^ (a), but with a vanishing edge 
region, may be constructed with pressure that is zero everywhere. 



VII. CONCLUSION 



We have developed a unified energy principle for analysing the MHD stability of plasmas consisting of both Taylor- 
relaxed and ideal MHD regions. The gauge a = ^xB for the vector potential, a, of linearized perturbations is used, 
with the equilibrium magnetic field B obeying a Beltrami equation, VxB /iB, in relaxed regions. A central result 
is that ^ obeys the same Euler-Lagrange equation whether ideal or relaxed MHD is used for perturbations, except 
in the neighbourhood of the magnetic surfaces where B-V is singular. To explore this difference, we developed a 
procedure for constructing global multi-region solutions in cylindrical geometry. Each region is a cylindrical shell in 
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which the plasma is either ideal MHD or Taylor-relaxed. The shells are nested and separated by ideal MHD barriers 
or interfaces of zero width. In this configuration, only Newcomb's small solutions are allowed for ideal MHD regions 
whereas in relaxed MHD only the odd-parity large solution and even-parity small solution are allowed. 

By way of illustration, we have evaluated stability implications using the different energy principles to a two interface 
plasma whose rotational transform profile resembles that of a reverse field pinch. A second motivation to study this 
configuration is the resolution of singular-limit problem encountered previously by Hole et al [\J ] in multi-region 
relaxed MHD, when two barrier interfaces become arbitrarily close. In that work, it was found that the stability of a 
two interface plasma with Taylor-relaxed regions does not approach the stability of a single interface configuration in 
the limit that the interface separation vanishes. 

We infer that the two interface configuration is unstable to a tearing mode, and that if the Taylor-relaxed region 
is replaced by an ideal MHD region of nonzero width, stability conclusions of the single interface and vanishing 
barrier region agree. By solving for the eigenvalues in the ideal treatment in the limit of vanishing barrier width, 
and comparing these to zero, we have put brackets on the pressure at field resonances. Using this, we have shown 
that providing the magnetic shear is monotonic across the plasma, it is always mathematically possible to construct 
a stable non-pathological pressure profile to connect a pressure jump. If the shear toggles sign however, such that 
an ideal resonance is in the plasma, stabilisation in the barrier region may require a pathological pressure profile (ie 
negative pressure). 

A conclusion of our work is that the stability of plasmas with jump in rotational transform across a barrier can 
not be concluded from a barrier treatment (such as Kaiser and Uecker or Hole et al) alone. One must also identify 
whether resonances exist in the barrier in the limit it is spatially resolved. If no resonances exist in the barrier and 
the plasma is stable, the structure of the pressure jump across the barrier in the limit it is resolved is irrelevant. If 
resonances exist in the barrier stability is then a function of the pressure profile within the barrier in the limit it is 
spatially resolved. 
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